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Abstract 

Contraction analysis is a stability theory for nonlinear systems where stability is 
defined incrementally between two arbitrary trajectories. It provides an alternative 
framework in which to study uncertain interconnections or systems with external in- 
puts, where it offers several significant advantages when compared with traditional 
Lyapunov analysis. Contraction-based methods are particularly useful for analyzing 
systems with uncertain parameters and for proving synchronization properties of non- 
linear oscillators. Existence of a contraction metric for a given system is a necessary 
and sufficient condition for global exponential convergence of system trajectories. For 
systems with polynomial or rational dynamics, the search for contraction metrics can 
be made fully algorithmic through the use of convex optimization and sum of squares 
(SOS) programming. The search process is made computationally tractable by relax- 
ing matrix definiteness constraints, whose feasibility indicate existence of a contraction 
metric, into SOS constraints on polynomial matrices. We illustrate the results through 
examples from the literature, emphasizing the advantages and contrasting the differ- 
ences between the contraction approach and traditional Lyapunov techniques. 

1 Introduction 

Contraction analysis is a stability theory for nonlinear systems where stability is defined 
incrementally between two arbitrary trajectories ^Hl- The existence of a contraction metric 
for a nonlinear system ensures that a suitably defined distance between nearby trajectories is 
always decreasing, and thus trajectories converge exponentially and globally. One important 
application of contraction theory is its use in studying the synchronization of nonlinear 
coupled oscillators j2S]- These oscillators present themselves in a variety of research fields 
such as mathematics, biology, neuroscience, electronics, and robotics. The use of coupled 
oscillators in each of these fields, as well as how contraction theory can be used to analyze 
networks of coupled identical nonlinear oscillators can be can be found in [23l and the 
references listed therein. 

Contraction theory nicely complements Lyapunov theory, a standard nonlinear stability 
analysis technique, as it provides an alternative framework in which to study convergence and 
robustness properties of nonlinear systems. For autonomous systems one can interpret the 
search for a contraction metric as the search for a Lyapunov function with a certain structure. 
This statement will be explained further in Section El There are, however, advantages to 
searching for a contraction metric instead of searching explicitly for a Lyapunov function. In 
particular, as we will show, contraction metrics are useful for analyzing uncertain nonlinear 
systems. In general, nonlinear systems with uncertain parameters can prove quite trou- 
blesome for standard Lyapunov methods, since the uncertainty can change the equilibrium 
point of the system in very complicated ways, thus forcing the use of parameter-dependent 
Lyapunov functions in order to prove stability for a range of the uncertain parameter values. 

Much of the literature on parameter-dependent Lyapunov functions focuses on linear 
systems with parametric uncertainty IH El C]- However, if a linear model is being used 
to study a nonlinear system around an equilibrium point, changing the equilibrium of the 
nonlinear system, necessitates relinearization around the new equilibrium. If the actual 
position of the equilibrium, in addition to the stability properties of the equilibrium, of 
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the nonlinear system depends on the uncertainty, it may be impossible to obtain any kind 
of closed form expression of the equilibrium in terms of the uncertain parameters. Thus, 
parameterizing the linearization in terms of the uncertainty may not be an option. 

A well-studied method of dealing with specific forms of nonlinearities is to model the 
nonlinear system as a linear system with bounded uncertainty. In particular, in |2] polytopic 
linear differential inclusions (LDIs), norm-bound LDIs, and diagonal norm-bound LDIs are 
considered. These techniques are computationally tractable as they reduce to convex opti- 
mization problems. Though these methods work for various kinds of uncertainty, it is also 
desirable to find methods to study the stability of nonlinear systems that do not easily admit 
linear approximations with the nonlinearities covered with uncertainty bounds. 

Contraction theory provides a framework in which to study the stability behavior of more 
general uncertain nonlinear systems. This framework eliminates many of the restrictions and 
problems that may be encountered when trying to analyze uncertain nonlinear systems with 
traditional linearization techniques or Lyapunov methods. This results from the fact that if a 
nominal system is contracting with respect to a certain contraction metric, it is often the case 
that the uncertain system with additive or multiplicative uncertainty within a certain range 
will still be contracting with respect to the same metric, even if the perturbation changes the 
position of the equilibrium of the system. Thus, it is possible to determine stability of the 
system for a range of values of the uncertain parameter without explicitly tracking how the 
uncertainty changes the location of the equilibrium. These ideas will be discussed further in 
Section 

Another interesting feature of the contraction framework is its relative flexibility in in- 
corporating inputs and outputs. For instance, to prove contraction of a class of systems with 
external inputs, it is sufficient to show the existence of a contraction metric with a certain 
structure. This feature, which will be discussed in Section |BJ is central in using contraction 
theory to prove synchronization of coupled nonlinear oscillators. 

To translate the theoretical discussion above into effective practical tools, it is desirable 
to have efficient computational methods to numerically obtain contraction metrics. Sum 
of squares (SOS) programming provides one such method. SOS programming is based on 
techniques that combine elements of computational algebra and convex optimization, and 
has been recently used to provide efficient convex relaxations for several computationally 
hard problems ^U]- In this paper we will show how SOS programming enables the search for 
contraction metrics for the class of nonlinear systems with polynomial dynamics. We discuss 
how to use SOS methods to find bounds on the maximum amount of uncertainty allowed in 
a system in order for the system to retain the property of being contracting with respect to 
the contraction metric of the unperturbed system. We also use SOS methods to optimize the 
contraction matrix search to obtain a metric that provides the largest symmetric uncertainty 
interval for which we can prove the system is contracting. 

This paper is organized as follows: in Section |21 we give background material on contrac- 
tion theory. Section El discusses sum of squares (SOS) polynomials and matrices. We present 
next an algorithm which uses SOS programming to computationally search for contraction 
metrics for nonlinear systems. We discuss why contraction theory is useful for studying 
systems with uncertain dynamics in Section El and external inputs in Section IHl Finally, in 
Section [71 we present our conclusions, and outline possible directions for future work. 
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2 Contraction Analysis 



Contraction analysis is a relatively recently developed stability theory for nonlinear systems 
analysis J^. The theory attempts to answer the question of whether the limiting behavior 
of a given dynamical system is independent of its initial conditions. More specifically, con- 
traction analysis is a theory in which stability is defined incrementally between two arbitrary 
trajectories. It is used to determine whether nearby trajectories converge to one another. 
This section summarizes the main elements of contraction analysis; a much more detailed 
account can be found in [T3]. 

We consider deterministic dynamical systems of the form 

i = f(x(t),t), (1) 

where f is a nonlinear vector field and x(t) is an n-dimensional state vector. For this analysis 
it is assumed that all quantities are real and smooth and thus that all required derivatives 
or partial derivatives exist and are continuous. This existence and continuity assumption 
clearly holds for polynomial vector fields. 

Under the assumption that all quantities are real and smooth, from equation ^ we can 
obtain the differential relation 

6Mt) = ^^{At),t)S^{t), (2) 

where (5x(t) is an infinitesimal displacement at a fixed time. For notational convenience from 
here on we will write x for x(t), but in all calculations it should be noted that x is a function 
of time. 

The infinitesimal squared distance between two trajectories is (5x^5x. Using the 
following equation for the rate of change of the squared distance between two trajectories is 
obtained: 

4(<^x^^x) = 25x^5x = 25x^-^5x. (3) 
dV ' (9x ^ ' 

If Ai(x, t) is the largest eigenvalue of the symmetric part of the Jacobian ^ (i.e. the largest 

eigenvalue of |(|^ + then it follows from Q that 



Integrating both sides gives 



^(^x^'^x) < 2Ai(x,t)5x'^5x. (4) 
at 

||5x|| < ||5xo||e^o^i(-.*)'^*. (5) 

If Ai(x, t) is uniformly strictly negative (i.e. (^ + ^) ^ V x, t), it follows from © 
that any infinitesimal length ||5x|| converges exponentially to zero. By path integration the 
distance of any finite path also converges exponentially to zero. 
A more general definition of length can be given by 

(5z^5z = 5x^M(x,t)(5x (6) 

where M(x, t) is a symmetric, uniformly positive definite and continuously different iable 
metric (formally, this defines a Riemannian manifold). This notion of infinitesimal distance 
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defined with respect to a metric can be use to define a finite distance measure between two 
trajectories with respect to this metric. Specifically, the distance between two points Pi and 
P2 with respect to the metric M(x, t) is defined as the shortest path length, in other words 
the smallest path integral Jp^ 5x^M(x, t)(5x, between these two points. Accordingly a ball 
of center c with radius R is defined as the set of all points whose distance to c with respect 
to M(x, t) is strictly less than R. 

Under the definition of infinitesimal length given in (jH)), the equation for its rate of change 
becomes 

^(5x^M5x) = 5x^(-!^ M + M-^ + M)5x (7) 
dt ax ox 

where M is shorthand notation for M(x, t). Convergence to a single trajectory occurs in 

regions where M + + M) is uniformly negative definite. It should be noted that 

M = M(x, t) = ^^Q^'^^ ^ + The above analysis leads to the following definition and 

theorem: 

Definition 1 ([13j). Given the system equations x = f(x, t), a region of the state space 
is called a contraction region with respect to a uniformly positive definite metric M(x, t) if 

M + + M) is uniformly negative definite in that region. 

Theorem 1 ( [13] ). Consider the system equations x = f(x, t). Assume a trajectory starts 
in a ball of constant radius that is defined with respect to the metric M(x, t), that is centered 
at a another given trajectory, and that is contained at all times in a contraction region with 
respect to the metric M(x, t). Then the first trajectory will remain in that ball and converge 
exponentially to the center trajectory. Furthermore, global exponential convergence to the 
center trajectory is guaranteed if the whole state space is a contraction region with respect to 
the metric M(x, t). 

Definition 1 provides sufficient conditions for a system to be contracting. Namely, the 
following should be satisfied: 

1. The matrix M(x, t) must be a uniformly positive definite matrix, i.e., 

M(x, t)yely0 Vx, t. (8) 

2. The metric variation + + M must be a uniformly negative definite matrix, 
i.e., 

R(x,t) = — M + M— + M ^ -el ^ Vx,t. (9) 
ox ox 

An explicit rate of convergence of trajectories f3 can be found by finding a M(x, t) that 
satisfies (jH)) and 

di'^ di ■ 

— M + M— + M ^ -/5M. (10) 
ox ox 

The notation above is standard; >~, and >z mean positive definite and positive semi definite 
respectively, while -< and ^ mean negative definite and negative semidefinite respectively. If 
the system dynamics are linear and M(x,t) is constant (i.e. M(x, t) = M), the conditions 
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above reduce to those in standard Lyapunov analysis techniques. Lyapunov theory shows 
that the system x(t) = A'x{t) is stable (i.e., all trajectories converge to 0) if and only if there 
exists a positive definite matrix M (i.e., M >- 0) such that A^M + MA -< 0. 

It should be noted that if a global contraction metric exists for an autonomous system, all 
trajectories converge to a unique equilibrium point, and we can always produce a Lyapunov 
function for the system from the contraction metric JH] ■ We assume, without loss of gener- 
ality, that the equilibrium is at the origin. If the system dynamics are f(x) and M(x) is a 
time-invariant contraction metric for the system, then ^(x) = f (x)-'"M(x)f (x) is a Lyapunov 

function for the system since y(x) > and y = f(x)^(g^M + + M)f(x) < -(3V. 
This shows that x = f(x) tends to exponentially, and thus that x tends towards a finite 
equilibrium point. 

For a constant metric M(x, t) = M, this reduces to Krasovskii's Method [H]. We note 
that for systems with uncertainty there are good reasons to search for a contraction metric 
to create Lyapunov function of this structure instead of searching for a Lyapunov function 
directly. These reasons will become clear in Section El 

The problem of searching for a contraction metric thus reduces to finding a matrix func- 
tion M(x, t) that satisfies the conditions above. As we will see, SOS methods will provide a 
computationally convenient approach to this problem. 

3 Sum of Squares (SOS) Polynomials and Programs 

The main computational difficulty of problems involving constraints such as the ones in (jHJ 
and © is the lack of efficient numerical methods that can effectively handle multivariate 
nonnegativity conditions. A convenient approach for this, originally introduced in ^7], is 
the use of sum of squares (SOS) relaxations as a suitable replacement for nonnegativity. We 
present below the basic elements of these techniques. 

A multivariate polynomial p{xi,X2, ■■■,Xn) = p(x) G ]R[x] is a sum of squares (SOS) if 
there exist polynomials /i(x), fm{'^) £ such that 

m 

Mx) = $^/^(x). (11) 

i=l 

The existence of a SOS representation for a given polynomial is a sufficient condition for its 
global nonnegativity, i.e., equation (fTTj) implies that ]9(x) > V x G M". The SOS condition 
(fTT|) can be shown to be equivalent to the existence of a positive semidefinite matrix Q such 
that 

p(x) = Z^(x)QZ(x) (12) 

where Z{x.) is a vector of monomials of degree less than or equal to deg(p)/2. This equivalence 
of descriptions between (fTTj) and (fT^ makes finding an SOS decomposition a computation- 
ally tractable procedure. Finding a symmetric positive semidefinite Q subject to the affine 
constraint (fT^ is a semidefinite programming problem ^71 ■ 

Using the notion of a SOS polynomial as a primitive, we can now introduce a convenient 
class of optimization problems. A sum of squares program is a convex optimization problem 
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of the form: 



J 



mm 



J 

subject to aifi + Cj 



is SOS for i 



1 



...,/ 



where the c/s are the scalar real decision variables, the wj^s are given real numbers that define 
the objective function, and the aij(x) are given multivariate polynomials. There has recently 
been much interest in SOS programming and SOS optimization as these techniques provide 
convex relaxations for various computationally hard optimization and control problems; see 
e.g. [m dl d 1201 and the volume ^. 

A SOS decomposition provides an explicit certificate of the nonnegativity of a scalar 
polynomial for all values of the indeterminates. In order to design an algorithmic procedure 
to search for contraction metrics, we need to introduce a similar idea to ensure that a 
polynomial matrix is positive definite for every value of the indeterminates. A natural 
definition is as follows: 

Definition 2 ([5J). Consider a symmetric matrix with polynomial entries S(x) G M[x]™'^'", 
and let y = [yi, . . . , ym]'^ be a vector of new indeterminates. Then S(x) is a sum of squares 
matrix if the scalar polynomial y"^S(x)y is a sum of squares in ]R[x, y]. 

For notational convenience, we also define a stricter notion: 

Definition 3. A matrix S(x) is strictly SOS z/S(x) — el is a SOS matrix for some e > 0. 

Thus, a strictly SOS matrix is a matrix with polynomial entries that is positive definite 
for every value of the indeterminates. An equivalent definition of an SOS matrix can be 
given in terms of the existence of a polynomial factorization: S(x) is a SOS matrix if and 
only if it can be decomposed as S(x) = T(x)-^T(x) where T(x) G ]R[x]^^™'. For example. 



is a SOS matrix for all values of a and k. Indeed, this follows from the decomposition 
M(x) = T(x)^T(x), where 



SOS matrices have also been used recently by Hoi and Scherer jS] and Kojima PU] to pro- 
duce relaxations of polynomial optimization problems with matrix positivity definiteness 
constraints. 

4 Computational Search for Contraction Metrics via 
SOS Programming 

As explained in Section |2l given a dynamical system, the conditions for a contraction metric 
to exist in regions of the state-space are given by a pair of matrix inequalities. In the case 
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of metrics M(x) that do not depend explicitly on time, relaxing the matrix definiteness 
conditions in © and © to SOS matrix based tests makes the search for contracting metrics 
a computationally tractable procedure. More specifically, the matrix definiteness constraints 
on M(x) (and R(x)) can be relaxed to SOS matrix constraints by changing the inequality 
M(x) — el ^ in (jSI) (where e is an arbitrarily small constant) to the weaker condition that 
M(x) be a strictly SOS matrix. With these manipulations we see that existence of SOS 
matrices M(x) and R(x) is a sufficient condition for contraction. 

Lemma 1. Existence of a strictly SOS matrix M(x) and a strictly SOS matrix — R(x) = 
— (^"^M + + M) is a sufficient condition for global contraction of an autonomous 

system x = f(x) with polynomial dynamics. 

Proof. By Theorem 1, a sufficient condition for contraction of any nonlinear system is the 
existence of uniformly positive definite M(x) and — R(x). A sufficient condition for uniform 
positive definiteness of M(x) and — R(x) is the existence of strictly SOS matrices M(x) and 
-R(x). □ 

This lemma can easily be extended to existence of certain SOS matrices implying con- 
traction with a convergence rate {3 by redefining R(x) as R(x) = g M + M£ + M + /3M. 
At this point, we do not know if the full converse of Lemma 1 holds. If a system is exponen- 
tially contracting, it is known that a contraction matrix always exists ^13j. Nevertheless, a 
system with polynomial dynamics may certainly be contracting under non-polynomial met- 
rics. Furthermore, even if a positive definite contraction matrix with polynomial entries M 
exists, it may not be the case that it is a SOS matrix. We notice, however, that some of 
these issues, such as the gap between "true" contracting metrics and SOS-based ones, can 
be bridged by using the more advanced techniques explained in |18j. 

4.1 Search Algorithm 

One main contribution of this work is to show how sum of squares (SOS) techniques can be 
used to algorithmically search for a time-invariant contraction metric for nonlinear systems 
with polynomial dynamics. Existence of a contraction metric for nonlinear systems certifies 
contraction (or convergence) of system trajectories. For systems with polynomial dynamics, 
we can obtain a computationally tractable search procedure by restricting ourselves to a 
large class of SOS-based metrics. 

As suggested by Lemma 1, the main idea is to relax the search for matrices that satisfy 
matrix definiteness constraints M(x) >- and — R(x) >- into SOS-matrix sufficient con- 
ditions. Equivalently, we want to find a polynomial matrix M(x) that satisfies SOS matrix 
constraints on M(x) and R(x). The SOS feasibility problem can then be formulated as 
finding M(x) and R(x) such that y"^M(x)y is SOS and — y"^R(x)y is SOS. 

More specifically, the detailed steps in the algorithmic search of contraction metrics for 
systems with polynomial dynamics are as follows: 

1. Choose the degree of the polynomials in the contraction metric, and write an affine 
parametrization of the symmetric matrices of that degree. For instance, if the degree is 
equal to two, the general form of M(x) is 
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aixf + a2XiX2 + asxl + 04X1 + 05X2 + ae bixf + 62X1X2 + b^xl + 64X1 + 65X2 + fte 
bixf + b2XiX2 + b^xl + 64x1 + 65x2 + fee cixf + C2X1X2 + 03X3 + C4X1 + C5X2 + ce 

where a,, ftj, and q are unknown coefficients. 

2. Calculate ^ and define R(x) := ^^M + M^ + M. Thus, R(x) will also be a symmetric 
matrix with entries that depend affinely on the same unknown coefficients a^, hi, and q. 

3. Change matrix constraints M(x) ^ Vx, and R(x) = H^M + + M ^ Vx into 
scalar constraints on quadratic functions p(x, y) = y^M(x)y > V x, y, and r(x, y) = 
y"^R(x)y = y-^(^"^M + + M)y < Vx, y, where y is an n x 1 vector of new indeter- 
minates. 

4. Impose SOS constraints onp(x, y), and — r(x, y), and solve the associated SOS feasibility 
problem. If a solution exists, the SOS solver will find values for the unknown coefficients, 
such that the constraints are satisfied. 

5. Use the obtained coefficients ai,bi,Ci to construct the contraction metric M(x) and the 
corresponding R(x). 

6. Optionally, for graphical presentation, independent verification, or if the convex optimiza- 
tion procedure runs into numerical error, further testing can be done to verify the validity 
of the computed solution. To do this, we can check if the matrix constraints M(x) >- 0, and 
R(x) -< hold over a range of the state space by finding and plotting the eigenvalues over 
this range. If a true feasible solution does not exist, the minimum eigenvalue of M(x) will 
be negative or the maximum eigenvalue of R(x) will be positive. Either one of these cases 
violates the matrix constraints which certify contraction. In most semidefinite programming 
solvers, the matrix Q in (|12p is computed with floating point arithmetic. If Q is near the 
boundary of the set of positive semidefinite matrices, it is possible for the sign of eigenvalues 
that are zero or close to zero to be computed incorrectly from numerical roundoff and for the 
semidefinite program solver to encounter numerical difficulties. Numerical issues are further 
discussed in Section 15.3.11 

7. An explicit lower bound on the rate of convergence can be found by using bisection to 
compute the largest /3 for which there exist matrices M(x) >- and R/3(x) = ^ M + M^ + 
M + pM^O. 

For the specific examples presented later in the paper, we have used SOSTOOLS, a 
SOS toolbox for MATLAB developed for the specification and solution of sums of squares 
programs [19j. The specific structure of SOS matrices, or equivalently, the bipartite form 
of the polynomials p(x, y) and r(x, y) is exploited through the option sparsemultipartite 
of the command sosineq that defines the SOS inequalities. Future versions of SOSTOOLS 
will allow for the direct specification of matrix SOS constraints. 
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We present next are two examples of using this procedure to search for contraction metrics 
for nonhnear systems with polynomial dynamics. The systems studied are a model of a jet 
engine with controller, and a Van der Pol oscillator. 



4.2 Example: Moore-Greitzer Jet Engine Model 

The algorithm described was tested on the following dynamics, corresponding to a Moore- 
Greitzer model of a jet engine, with stabilizing feedback operating in the no-stall mode pjj. 
In this model, the origin is translated to a desired no-stall equilibrium. The state variables 
correspond to0 = $ — l,'?/' = ^~ ^co ~ 2, where $ is the mass flow, \E' is the pressure rise 
and \l/co is a constant JT]. The dynamic equations take the form: 



' " 









3(f) 



(13) 



The only real-valued equilibrium of the system is (p = 0, ip = 0. This equilibrium is stable. 

The results of the algorithmic search for SOS matrices M(x) and — R(x) of various 
orders are given in Tabled Values in the table, except the final row, are output values from 
SeDuMi j22j, the semidefinite program solver used as the optimization engine in solving 
the SOS program. CPU time is the number of seconds it took for SeDuMi's interior point 
algorithm to find a solution. As expected, the computation time increases with the degree 
of the polynomial entries of M(x). Feasibility ratio is the final value of the feasibility 
indicator. This indicator converges to 1 for problems with a complementary solution, and to 
— 1 for strongly infeasible problems. If the feasibility ratio is somewhere in between, this is 
usually an indication of numerical problems. The values pinf and dinf detect the feasibility 
of the problem. If pinf = 1, then the primal problem is infeasible. If dinf = 1, the dual 
problem is infeasible. If numerr is positive, the optimization algorithm (i.e., the semidefinite 
program solver) terminated without achieving the desired accuracy. The value numerr = 1 
gives a warning of numerical problems, while numerr = 2 indicates a complete failure due to 
numerical problems. 

As shown in Table ^ for this system no contraction metric with polynomial entries of 
degree or 2 could be found. This can be certified from the solution of the dual optimization 
problem. Since SeDuMi is a primal-dual solver, this infeasibility certificates are computed 
as a byproduct of the search for contraction metrics. 

An explicit lower bound for the rate of convergence of the trajectories of the jet engine 
model, i.e., the largest value P for which matrices M(x) >- and R/3(x) = ^ M + + 
M + /3M ^ were found, was f] = 0.818. 

We remark that for this system, it is also possible to prove stability using standard Lya- 
punov analysis techniques. However, we illustrate stability of this example from a contraction 
viewpoint because contraction theory offers a good approach to study this system when there 
is parametric uncertainty in the plant dynamics or feedback equations. For example, in the 
no-stall mode, the jet dynamics equations are 













—u 



(14) 
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Degree of polynomials in M(x) 





2 


4 


6 


CPU time (sec) 


0.140 


0.230 


0.481 


0.671 


Feasibility ratio 


-1.000 


-0.979 


1.003 


0.990 


pinf 


1 


1 








dinf 














numerr 





1 








M ;^ 0, i? ^ conditions met? 


no 


no 


yes 


yes 



Table 1: Contraction matrix search results for closed-loop jet engine dynamics. 

where m is a control variable. If a nominal stabilizing feedback control u can be found 
(e.g., using backstepping [TT] or some other design method), the SOS techniques described 
in Section 15.11 provide a way to find other stabilizing feedback controls which are centered 
around the nominal control. For example, if a stabilizing linear feedback control u = kicj) + 
k2il) can be found, we can interpret ki and k2 as uncertain parameters and use the methods 
described in Section 15.11 to search for ranges of gain values centered around the nominal 
values ki and ^2 that will also stabilize the system. 

4.3 Example: Van der Pol Oscillator 

A classic example that has played a central role in the development of nonlinear dynamics 
is given by the Van der Pol equation 

X + a{x^ + k)x + uPx = 0, (15) 

with a > 0, A;, and u as parameters. Historically this equation arose from studying nonlinear 
electric circuits used in the first radios When k < 0, the solutions of (fT3j) behave 

like a harmonic oscillator with a nonlinear damping term a(x^ + k)x. The term provides 
positive damping when |a;| > A; and negative damping when \x\ < k. Thus, large amplitude 
oscillations will decay, but if they become too small they will grow larger again [21] . li k > 
all trajectories converge to the origin. 

In Table 121 we present the results of running the contraction matrix search algorithm for 
the system 



Xi 




X2 


X2 




—a{xl + k)x2 — oJ^X\ 



with a = 1,0; = 1, which is the state-space version of the Van der Pol oscillator (|15|). We 
present solution for various values of A;, with a contraction matrix with entries that are 
quartic polynomials. 

As a natural first step we searched for a constant contraction metric. None could be 
found algorithmically. This was expected as it is easily shown analytically that a constant 
contraction matrix for this system does not exist. If M is constant, then 



M = 



a b 


di 


1 


b c 


' 5x 


—1 — 2xiX2 —xf — k 



—2b — AbxiX2 a — bx\ — kb — c — 2cx\X2 

a — bxi — kb — c — 2cxiX2 2b — 2cx\ — 2kc 
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Degree of polynomials in M(x) 


4 


4 


4 


4 


4 


4 


4 


4 


4 


4 


k 


-10 


-1 


-0.1 


-0.01 


-0.001 


0.001 


0.01 


0.1 


1 


10 


pinf 


1 


1 


1 























dinf 
































numerr 


1 


1 


1 


1 


1 

















M y 0, R < conditions met? 


no 


no 


no 


no 


no 


yes 


yes 


yes 


yes 


yes 



Table 2: Contraction matrix search results for oscillator dynamics. 




(a) k = 0.5 (b) k = -0.5 



Figure 1: Phase plots of Van der pol Oscillator. 

For R to be negative definite Rn must be negative for all values of Xi, X2- In other words 
—2b — 46x1X2 < or —1 < 2x1X2. This clearly does not hold for all values of Xi, and X2. A 
more complicated analysis (or a duality argument) also shows why there is no contraction 
matrix with quadratic entries for this system. 

The algorithm finds a contraction function for the system x -|- (x^ -|- k)x + x = when 
k > but not when k < 0. As shown in Figure ^ the trajectories of the oscillator converge 
to zero when k > 0, and converge to a limit-cycle when k < 0. Thus, the results of the 
contraction metric search is as expected. Since all trajectories converge to the origin when 
> we expect that a contraction metric exists for the system. In the case where A; < 
the origin is an unstable fixed point and thus the system is not contracting. 

Since for k < the system is not contracting, we should not be able to find a contraction 
function. It should be noted that the converse does not hold. The fact that we cannot find 
a contraction function does not necessarily mean that the system is not contracting. This is 
because finding an SOS representation of the constrained quadratic functions is a sufficient 
condition for their positivity, not a necessary one. 

It should be noted that for the example above, we can prove stability through Lyapunov 
analysis, and SOS programming can also be used to find Lyapunov functions [ISj. However, 
we illustrate this example here as contraction theory applied to a slightly modified version of 
this system provides a nice way to prove synchronization of coupled Van der Pol oscillators. 
This will be discussed in Section El This synchronization property is much more difficult to 
prove with standard Lyapunov methods. 
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5 Contraction Metrics and Systems with Uncertain 
Dynamics 

5.1 Uncertainty Analysis with Contraction Metrics and SOS Pro- 
gramming 

From the robust control perspective, one of the most appeahng features of contraction the- 
ory is the fact that it provides a natural framework in which to study uncertain nonlinear 
systems where the parametric uncertainty changes the location of the equilibrium points. In 
general, standard Lyapunov analysis does not handle this situation particularly well, since 
the Lyapunov function must track the changes in the location of the steady-state solutions, 
thus forcing the use of parameter-dependent Lyapunov functions. However, in general it may 
be impossible to obtain any kind of closed form expression of the equilibria in terms of the 
parameters, thus complicating the direct parametrization of possible Lyapunov functions. 

Much attention has been given to robust stability analysis of linear systems (e.g., [HHHEl 
El El])- Less attention, however, has been paid to nonlinear systems with moving equilibria. 
Two papers addressing this issue are [Elllj. The approach in ^14J is to consider systems 
described by the equations 



where x is a real n-vector, f and h are continuously differentiable functions, and h(x) 
represents the uncertainties or perturbation terms. Given an exponentially stable equilibrium 
Xe, il4|| establishes sufficient conditions by using the linearization of the system to produce 
Lyapunov functions which prove existence and local exponential stability of an equilibrium 
Xe for (fT7|) with the property |xe — Xe| < e where e is sufficiently small. 

Since the approach in is essentially based on a fixed Lyapunov function, it is more 
limited than our approach using contraction theory and SOS programming, and can prove 
stability only under quite conservative ranges of allowable uncertainty. By "allowable" we 
mean that if the uncertainty is in this range, the equilibrium remains exponentially stable 
under the uncertainty. A quantitative measure of this conservativeness will be given in 
Section 15.2.11 where we discuss the results of the method apphed to an uncertain model of 
a Moore-Greitzer jet engine and compare them to an approach via contraction theory and 
SOS programming. 

The approach in ^ is to linearize the dynamics around an equilibrium which is a function 
of the uncertain parameter (xq = g(5), 6 E Q) and then use structured singular values to 
determine the eigenvalues of the linearized system ^ = A(5)z if A(S) is rational in 5. If 
A{6) is marginally stable, no conclusions can be made about the stability of the nonlinear 
system. 

The contraction theory framework eliminates the need for linearization, and even the 
need to know the exact position of the equilibrium, in order to analyze stability robustness 
in uncertain nonlinear systems. In contrast to the Lyapunov situation, when certain classes 
of parametric uncertainty are added to the system, a contraction metric for the nominal 
system will often remain a contraction metric for the system with uncertainty, even if the 
perturbation has changed the equilibrium of the nonlinear system. 
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As noted in Section |2l if a global time-invariant contraction metric exists for an au- 
tonomous system, and an equilibrium point exists for the system, all trajectories converge 
to a unique equilibrium point, and we can always produce a Lyapunov function of the form 
y(x) = f (x)"^M(x)f (x). When a system contains parametric uncertainty, this formula yields 
the parameter-dependent Lyapunov function V^(x, S) = f (x, (5)"^M(x)f (x, 6) for ranges of the 
parametric uncertainty 6 where the contraction metric for the nominal system is still a con- 
traction metric for the system with perturbed dynamics. Thiis, if a contraction metric can 
be found for the system, we can easily construct a Lyapunov function which tracks the 
uncertainty for a certain range. 

5.1.1 Case 1: Bounds on the uncertainty range for which the system remains 
contractive with respect to the nominal metric. 

We can estimate the range of uncertainty under which the contraction metric for the nominal 
system is still a contraction metric for the perturbed system. To calculate this range, a SOS 
program can be written to minimize or maximize the amount of uncertainty allowed subject 
to the constraint Rsix) = ^'^M + M|| + M(f5(x)) -< 0, where f5(x) are the dynamics for 
the system with parametric uncertainty. The uncertainty bound is a decision variable in the 
SOS program and enters the constraint above in the ^ and isi'x.) terms. 

If we have more than one uncertain parameter in the system, we can find a polytopic inner 
approximation of the set of allowable uncertainties with SOS Programming. For example, if 
we have two uncertain parameters, we can algorithmically find a polytope in parameter space 
for which the original metric is still a contraction metric. The convex hull of four points, 
each which can be found by entering one of the four combinations, {Si, 62) = (7, 7), ^2) = 
(7, -7), (61,62) = (-7,7), or {61,62) = (-7, -7), into the uncertainty values in fs=[Su52]T i^) 

and then maximizing 7 subject to the constraint R^(x) = ^ M + + M(f^(x)) ^ 0, 

defines a polytope over which stability is guaranteed. 

5.1.2 Case 2: Search for a contraction metric that guarantees the largest sym- 
metric uncertainty interval for which the system is contractive. 

Alternatively, we can instead optimize the search for a metric M(x) that provides the largest 
symmetric uncertainty interval for which we can prove the system is contracting. If the scalar 
uncertainty 6 enters the system dynamics affinely, in other words if f(x) = fi(x) + 5f2(x), 
we can perform this optimization as follows. First write R(x, 5) = Ro(x) + (5Ri(x). To 
find the largest interval (—7, 7) such that for all 6 that satisfy —7 < 5 < 7 the system is 
contracting, introduce the following constraints into an SOS program: 

M(x) y 0, Ro(x) + 7Ri(x) -< 0, Ro(x) - 7Ri(x) -< 0. 

We note that 7 multiplies the scalar decision coefficients Oj, bi, and Cj in Ri(x) and thus we 
must use a bisection procedure to find the maximum value of 7 for which there exists SOS 
matrices M(x), Ro(x) and Ri(x) that satisfy the constraints above. 

If there are two uncertain parameters that enter the system dynamics affinely, we can 
extend the procedure above as follows: To find the largest uncertainty square with width 
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Degree of polynomials in M(x) 


4 


6 


6 range 


(-0.126,0.630) 


( -0.070, 0.635) 



Table 3: Range of perturbation where closed-loop uncertain jet engine dynamics given in 
(|19p are contracting with respect to the nominal metric. 



and height 7 such that for all and 62 that satisfy —7 < 5i < 7 and —7 < ^2 < 7 the 
system is contracting, first write R(x, 61, 62) = Ro(x) + 5iRi(x) + 52R'2(x), Then introduce 
the following constraints into and SOS program: 

M(x) y 0, Ro(x) + 7Ri(x) + 7R2(x) -< 0, Ro(x) + 7Ri(x) - 7R2(x) -< 

Ro(x) - 7Ri(x) + 7R2(x) ^ 0, Ro(x) - 7Ri(x) - 7R2(x) ^ 0. (18) 

Next, as in the scalar uncertainty case, use a bisection procedure to find the maximum value 
of 7 for which there exists SOS matrices M(x), Ro(x), Ri(x) and R2(x) that satisfy the 
constraints above. In the case of a large number of uncertain parameters, standard relaxation 
and robust control techniques can be used to avoid an exponential number of constraints. 



5.2 Example: Moore- Greitzer Jet Engine Model with Uncertainty 
5.2.1 Scalar Additive Uncertainty 

As described above, SOS programming can be used to find ranges of uncertainty under which 
a system with uncertain perturbations is still contracting with the original contraction metric. 
The contraction metric found for the deterministic system continues to be a metric for the 
perturbed system over a range of uncertainty even if the uncertainty shifts the equilibrium 
point and trajectories of the system. For the Moore-Greitzer jet engine model, the dynamics 
in (fT^ were perturbed by adding a constant term 6 to the first equation. 



" " 









+ 6 



(19) 



In Table |3] we display the ranges of d where the system was still contracting with the 
original contraction metric for 4*^^ and 6^^ degree contraction metrics. Note the range of 
allowable uncertainty is not symmetric. 

When instead we optimized the contraction metric search to get the largest symmetric S 
interval we obtained the results listed in Tabled A 6*^ degree contraction function finds the 
uncertainty range \6\ < 1.023. Because a Hopf bifurcation occurs in this system at 5 ~ 1.023, 
making the system unstable for 6 > 1.023, we can conclude that the 6th degree contraction 
metric is the highest degree necessary to find the maximum range of uncertainty for which 
the system is contracting. The Hopf bifurcation is shown in Figure |21 

Using the techniques in [T^ we computed the allowable uncertainty range for the system 
given in (jl9|) as |5| < 5.1 x 10~^. In the notation of ^3], we calculated the other parameters 
in Assumption 1 of [Hj as: h = [S, 0]^, \A~^\oo = 1, |-Dh(xe)|oo = 0) = 
|h(xe)|oo = where 6 is the perturbation term in (fT^ . The allowable range \6\ < 1.023 
computed via contraction theory and SOS programming is much larger than the allowable 
uncertainty range \5\ < 5.1 x 10~^ computed with the techniques in [Hj. 
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Degree of polynomials in M(x) 


4 


6 


8 


S range 


\S\ < 0.938 


\S\ < 1.023 


\S\ < 1.023 




Figure 2: Hopf bifurcation in uncertain jet dynamics. 



5.2.2 Scalar Multiplicative Uncertainty 

The approaches in Section also apply to multiplicative uncertainty, since the multiplica- 
tive coefficients enter affinely in the constraints in the SOS program. Tables El and IHl present 
the results of the described uncertainty analysis on the following system, which is equation 
(jl3|) with multiplicative uncertainty. 



" " 













(20) 



5.2.3 Multiple Uncertainties 

We consider next the system that results from introducing two additive uncertainties to the 
jet dynamics in equation (fT^ . We computed an uncertainty polytope (shown in Figure 
for which the system 

J 

is guaranteed to be contracting with respect to the original metric. Table [7| shows the results 
of optimizing the contraction metric to find the largest uncertainty square with width and 
height 7 such that for all 6i and S2 that satisfy —7 < (5i < 7 and —7 < 52 < 7 the system is 
contracting. 



Degree of polynomials in M(x) 


4 


6 


6 range 


(0.9767, 5.8686) 


(0.9796, 3.9738) 



102 _ 103 + 

3(p- ip + 52 



(21) 



Table 5: Range of perturbation for which the uncertain system given in ()2()|1 is contracting 
with respect to the nominal metric. 
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Degree of polynomials in M(x) 


4 


6 


8 


S range 


(1 -0.247,1 + 0.247) 


(1 - 0.356,1 + 0.356) 


(1 - 0.364, 1 + 0.364) 



Table 6: Symmetric range of perturbation where uncertain closed-loop jet engine dynamics 
given in ()2U|) are contracting. 




-0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 

52 



Figure 3: Polytopic region of uncertainty where closed-loop jet engine dynamics given in 
(|2ip are contracting with respect to nominal metric. 



Degree of polynomials in M(x) 


4 


6 


8 


Height and width of allowed uncertainty box 


0.7093 


0.7321 


0.7346 



Table 7: Symmetric range of perturbation where uncertain closed-loop dynamics given in 
(PT|) are contracting. 
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6 Contraction Metrics for Systems with External In- 
puts 



6.1 Stability Analysis of Systems with External Inputs 

Another interesting feature of the contraction framework is the relative flexibihty in incor- 
porating inputs and outputs. For instance, to prove contraction of a class of systems with 
external inputs, it is sufficient to show the existence of a polynomial contraction metric with 
a certain structure. This is described in the following theorem. 



Theorem 2. Let 

Xi = fi(xi, X2, ■ ■ ■ , Xn) 

Xk = fk{Xl,X2,. . . ,Xn) 

Xk+1 = fk+l{xi,X2,...,Xn) +Vk+l{u) 



Xn = fniXi,X2,...,Xn)+Vn{u) (22) 

be a set of nonlinear coupled differential equations where only the last n — k depend explicitly 

on u{t). If there exists anxn matrix M(xi, . . . , Xk) such that M and M+|^ M+M^ -< 
then the system is contracting for all possible choices ofu{t). 

Proof. For not at ional convenience, let xi = [ Xi ... Xk = fi(xi,X2) andx2 = [ Xk+i ... 
f2(xi,X2,M). The metric M(a:i, 0:2..., x^) = M(xi) is independent of X2, and thus = 
V2,i. Since ^ also vanishes, it follows that ^i,j, M,^ = + + ff = 

^dxl^ ^it • Thus M(xi) is not a function of u{t). In addition, ^ has no dependence on u{t) 
because f (x, u) = h(x) + v(m). Thus, if there exists anxn matrix M(xi, Xk) such that 
M >- and M + ^ M + M^-<0, then the system in (j^^ is contracting for any value of 

U{t). XX ^ 

Through the example considered in the following section, we will illustrate how TheoremEl 
is particularly useful in proving synchronization of nonlinear oscillators, an issue explored in 
more detail in [22] • Theorem |21 can be easily extended to the case where u{t) is a vector (i.e. 

U{t) = [ui{t),...,Um{t)f). 



6.2 Coupled Oscillators 

Contraction Theory is a useful tool to study synchronization behaviors of various config- 
urations of coupled oscillators. For simplicity, we only consider here the case a pair of 
unidirectionally coupled oscillators; more complicated and general couplings are discussed in 
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A state-space model of two unidirectionally coupled oscillators (only one oscillator influ- 
ences the other) is 

± = f(x,t) 

y = f(y,^) + u(x) -u(y), 



(23) 



where x, y G M™", are the state vectors, f (x, t) and f (y, t) are the dynamics of the uncoupled 
oscillators, and u(x) — u(y) is the coupling force^. The following theorem is a slightly 
modified version of Theorem 2 in 



Theorem 3. If y = f(y) + u(y) — u(x) in h2!^) is contracting with respect to y over the 
entire state space for arbitrary u(x)^, the two systems will reach synchrony (i.e. y{t) and 
x(t) will tend toward the same trajectory) regardless of initial conditions. 

Proof. The system y = f (y) — u(y) + u(x) with input u(x) is contracting with respect to y 
over the entire state space and y{t) = x(t) is a particular solution. Thus, by the properties 
of contraction, all solutions converge exponentially to y{t) = x(t). □ 

Theorem |21 becomes especially powerful when the vector field appearing in the second 
subsystem of has the structure described in equation (j22I)^- We illustrate this in the 
next example. 



6.3 Example: Coupled Van der Pol Oscillators 

Consider two identical Van der Pol oscillators coupled as 



X + q;(x^ + k)x + uP'x 
ij + a{y^ + k)y + u^y 





ari{x — y) 



(24) 



where a > 0, a; > 0, A; are arbitrary constants. We note that if A; < 0, trajectories of the 
individual oscillator dynamics converge to a limit cycle. See Figure |l(b)| We first write 
these coupled systems in state-space form to get the equations in the form of (j23p . Their 
state-space form is 











Xi 




< 



















X2 



—a{xl + k)x2 — uP'Xx 



Z/2 

-a{y\ + A; + 77)2/2 - ^'^Vx + 01^X2 



(25) 



where a > 0, w > 0, fc are arbitrary constants. 



^An example of coupled oscillators whose state-space representation is in this form is 

{X + aix^ + k)x ijj^ X =0 
y + a{y^ + k)y + iv^y ^ ar]{x - y) 
^By contracting with respect to y for arbitrary u(x) we mean that the system y = f(y) — u(y) + u(x), 
where y is the state vector and u(x) is an arbitrary driving function, is contracting for all inputs u(x). 

■^If it does not have such a structure and u(x) drives each component of y, we lose degrees of freedom in 
the possible forms of our contraction metric. If u(x) drives each component of y the only possible contraction 
metric is a constant. 
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By Theorem |2l this pair of unidirectional oscillators will reach synchrony regardless of initial 
conditions if 



y = f(y) -u(y) + u(x) 



?/2 

-a{yj + k + r])y2 - w'^yi + ar]X2 



(26) 



is contracting with respect to y for arbitrary values of u(x) = X2- We see by Theorem |21 
that for this to occur, we must find a contraction metric M(y) that is only a function of yi 
(i.e. M(y) = M(|/i)). 

When the search algorithm described in Section 14.11 was applied to find a metric that 
satisfied M(y) = M(?/i) as well as M(y) >- and R(y) -< 0, none were found. However, it 
is shown in the appendix, which is a modified version of the appendix of [23j, that a metric 
that satisfies M(y) >- and R(y) ^ implies asymptotic convergence of trajectories of 
system fl2fij) . A system with this metric that satisfies M(y) >- and R(y) ^ is called 
semi- contracting [iniESI- 
The metric 



M(y) 



uj"^ + a^iyl + k + r]f 
a{yl + k + r]) 



a{yl + k + r]) 
1 



(27) 



that appears in [23 is only a function of yi and satisfies M(y) >- and R(y) ^ for the 
system dynamics (pUj) if a > and {k + rj) > 0. For this M and the system equation 
we have 



R=M+— M+M— 
oy oy 



-2auj'^yl - 2auj'^{k + 77) 









(28) 



For a > 0, (A; + 77) > 0, M(y) ^ and R(y) ^ 0. Since (j2I|) and (j2Hl) show analytically that 
the system ()2(ij) is semi-contracting we used our search algorithm to search for a metric with 
M(y) ^ and R(y) ^ 0. 



6.3.1 Search for a Semidefinite R Matrix: Numerical Problems and Solutions 

A minor problem that one may encounter when searching for contraction metrics, depending 
on the structure of polynomial constraints, is that the resulting optimization problem may 
be feasible, but not strictly feasible. This can cause numerical difficulties in the algorithms 
used in the solution procedure. In many cases, however, this can be remedied by a intro- 
ducing a presolving stage in which redundant variables are eliminated. When we ran the 
search algorithm based on Theorem |21 and only searched for M as a function of yi, no valid 
solution was found even if we only constrained R to be negative semidefinite and not strictly 
negative definite. Since the analytic solution (j28|) was feasible but not strictly feasible, we 
hypothesized there was numerical error in the algorithm. Based on knowledge of the analytic 
solution (|^. we thus constrained R22 = and R12 = 0, eliminated redundant variables, and 
then searched for a solution in the resulting lower dimensional space^. With these constraints 
in place, a solution was found with the search algorithm. 

''Setting R22 = 0, and R12 — leads to redundant decision coefficients in tlie polynomial entries of M 
and R. If these redundant variables are eliminated through a presolving stage, the search algorithm finds 
M and R ^ 0. 
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7 Conclusions 

In this paper we have described how SOS programming enables an algorithmic search for 
contraction metrics for the class of nonlinear systems with polynomial dynamics. We also 
have illustrated the results through several examples. 

These examples illustrate how contraction analysis offers several significant advantages 
when compared with traditional Lyapunov analysis. Contraction analysis provides relative 
flexibility in incorporating inputs and outputs. It is also particularly useful in the analysis of 
nonlinear systems with uncertain parameters where the uncertainty changes the equilibrium 
points of the system. It is often the case that if the nominal system is contracting with 
respect to a metric, the uncertain system with additive or multiplicative uncertainty will 
still be contracting with respect to the original metric, even if the perturbation changes the 
equilibrium of the system. In addition, a slightly modified version of the standard algorithmic 
search allows us to optimize the search to obtain a contraction metric that provides the largest 
uncertainty interval for which we can prove the system is contracting. 

Subjects of future research include a careful evaluation of how the computational re- 
sources needed by the algorithm scale with system size, as well as the benefits and limitations 
of this approach in the context of other nonlinear system analysis techniques. 



A Proving Asymptotic Convergence of Coupled Van 
der Pol Oscillators With a Negative Semidefinite R 
Matrix. 



R(y) 



This appendix is a modified version of the appendix in j23|. Consider the system given in 
(I26|) . Consider a 2 x 2 matrix M(y) that is uniformly positive definite, and a corresponding 
R(y) matrix that is uniformly negative semidefinite, but not uniformly negative definite. 
Since is a two-dimensional system, we can assume without loss of generality that R(y) 
is of the form 

• -K{y) 


where K{y) > V y. Let <5y = ( 6yi 6y2 )^ = { Sy 6y )^ . where yi and y2 are the 
variables in equation ()26|) . With this R(y) and M(y) matrices, the general definition of 
differential length given in (jH)) and associated equation for rate of change of length ((Tj) are 

5z^5z = 5y^M(y)5y 

and 

|.z-.z ^ |(.y-M(y).y) 
= 5y^R(y)5y 

= -K{y)6yl (29) 
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Since J^(5y^M(y)5y) < and b-L^b^ > 0, Sz^Sz has a limit at t goes to infinity. We 
will prove through a Taylor series argument that all trajectories of this system converge 
asymptotically. If 6y = 6yi ^ 0, then 

Sz^dz{t + dt) - 6z^6z{t) = -K{y){6yifdt + 0{{dtf) 

while if 5yi = 0, 

dz^dz{t + dt) - 5z^5z(t) = -2K{y){dy2f^ + 0{{dtf). 

Since 5z^5z converges, 5z^ 5z{t+ dt) — 5z^ 5z{t) approaches zero asymptotically and hence 5yi 
and 5y2 or equivalently 5y and 5y both tend to zero. Thus, for any input u(x) all solutions of 
system fl2(j|l converge asymptotically to a single trajectory independent of initial conditions, 
and the unidirectional oscillators given in will reach synchrony asymptotically regardless 
of initial conditions. 
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